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Abstract 

We review the current status of the global DGLAP analysis of nuclear parton distribution functions, nPDFs, fo- 
cusing on the recent EPS09 analysis (T), whose output, EPS09NLO, is the best-constrained NLO nPDF set on the 
market. Collinear factorization is found to work very well in the kinematical region studied. With the error sets 
released in the EPS09 package one can compute how the nPDF-related uncertainties propagate into factorizable nu- 
clear hard-process cross sections. A comparison with the other existing NLO nPDF sets is shown, and the BRAHMS 
forward-77 hadron data from d+Au collisions are discussed in the light of the EPS09 nPDFs and their error sets. 

1 Introduction 

In this talk, I will discuss the nPDFs which are obtained through genuinely global analyses, which are based on the 
DGLAP evolution (2) and where one uses various types of nuclear hard process data in extracting the nPDFs. Such 
analyses are a test of perturbative QCD (pQCD) and collinear factorization, schematically expressed as 

cT A B^h+x = ]Ti^(Q 2 )®/f {Q 2 )®^ h+x ,{Q 2 ) + 0{l/Q 2 ), (1) 

where the factorization/renormalization scales are sufficiently large, Q 2 ^> Aq CD , and the PDFs f^, f? are universal 
and transferable from one hard process to another. 

Traditionally, both for the free and bound proton PDFs, the global DGLAP analysis procedure has been as follows: 
one starts with the non-perturbative input, an initial assumption, for the PDFs fi(x, Qq, {a,}), which are expressed 
in terms of a set of fit parameters {a;} and where the initial scale Qq typically is of the order of 1 GeV. QCD sum 
rules are imposed at this stage. The PDFs at larger scales, fi(x, Q 2 > Qq, {cti}), are obtained through the standard 
pQCD DGLAP evolution to the required order (LO, NLO,...). After this, one compares the various hard-process cross 
sections, which are computed using the set {fi(x, Q 2 , {ai})}, with existing data. The best set of PDFs is determined 
through iteration in the a^-parameter space, by finding the set {a^} which minimizes the \ 2 of the fit. After this, one 
has to perform a rather non-trivial error analysis to determine the non-correlated "error PDF sets" which allow one to 
study the propagation of the nPDF-uncertainties into hard cross sections. Thus, the outcome of such global analysis is 
the best PDF set {fi(x,Q 2 > Q 2 ,)}, supplemented with a number of error sets. 

Here I will not be addressing the QCD origin of the nuclear effects in PDFs - for discussion on this important 
question see e.g. Refs. J3] |4j [5] |6) . I will focus on the linear DGLAP evolution only, studies of other evolution 
equations (BFKL, BK,...), saturation or nonlinear effects [7| will not be addressed here either. As indicated by Eq. ([lj, 
power corrections (see e.g. Refs. [8 9|) will be ignored also. 

The global analyses of PDFs and nPDFs are quite challenging, since in general the data lie in correlated regions of 
x and Q 2 , and since there are statistical, systematic and also additional normalization errors to deal with but no unique 
way of propagating these uncertainties into the PDFs. As the parameter space is 15-30 dimensional, very fast DGLAP 
and cross-section solvers are needed especially in the NLO (for technical details, see ifTUl ). 

The global DGLAP analyses have resulted in excellent fits for collinearly factorized free proton PDFs: sets like 
CT10 03], MSTW H2l, and (from neural networks) NNPDF2.0 OS] are nowadays available. The global analysis 
of nPDFs, however, is even more challenging than that of the free proton PDFs, since in addition to the x and Q 2 
dependence one needs to consider also the A dependence of the PDFs. There also is a smaller amount of data from 
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fewer types of processes available as constraints, and the kinematical coverage of the nuclear hard-process data is 
more limited than in the free proton PDF case. 

Genuinely global analyses of nPDFs have been performed now over a decade. The pioneering LO analysis, where 
the x 2 minimization was made by eye, resulted in the nPDF set EKS98 04), which still is valid and useful. The first 
error estimates on nPDFs were provided for the LO sets HKM ifTBI and HKN04 [16|. The first NLO nPDFs were 
given in the nDS set [17| and the first error estimates in the NLO case in the HKN07 set |18]. The latest progress was 
made in the EPS09 analysis HI which includes RHIC data and now for the first time provides, in addition to the best 
NLO fit, 30 concrete nPDF error sets (for each A). With EPS09, which currently represent the state-of-the-art NLO 
nPDFs, one has reached the same degree of technical rigor in the global analysis for nPDFs as achieved already some 
time ago in the free proton case. 



2 EPS09 global analysis framework 

Let me next describe the main features in the EPS09 analysis framework, for details, see the Ref. [1 1. We define the 
bound proton PDFs relative to the free proton ones, taking the set CTEQ6.1M [ 19] as our baseline: 

tf(z, Q 2 ) EE R?(X, Q 2 )/fTEQ6.1M (X; g2); (2) 

where i is the parton flavour. Choosing CTEQ6.1M means that we work in the A1S and zero-mass variable flavor- 
number schemes. Isospin symmetry is assumed, so that the average u-quark distribution in a nucleus A becomes 
ua{x, Q 2 ) = ^Ju( x j Q 2 ) + A ~s L Z fd( x ' Q 2 )> an( l similarly f° r the <i-quark and u, d antiquarks. 

As illustrated in Fig. 1 (left), one has to deal with the following x-, Q 2 - and A-dependent nuclear effects in the 
nPDFs: the "shadowing" depletion observed in the DIS structure function ratios \f£ j \F® at x <C 0.1, the "anti- 
shadowing" excess around x ~ 0.1, the "EMC-effect" depletion at 0.2 < x < 0.7 and the "Fermi motion" excess 
towards x = 1 and beyond. These effects are in the nPDFs at our (CTEQ6.1M's) initial scale Qq — 1.69 GeV 2 . In 
the lack of sufficient data constraints, we are forced to start with only three different modification ratios: Rq(x, Qq), 
Ry(x, Qq), and Rs{x, Qq) for gluons, valence quarks and sea quarks, correspondingly. The A dependence is em- 
bedded in the A dependence of the parameters, such as the one controlling the antishadowing: y a = y a (C)(A/12) Pa , 
with carbon as the reference nucleus. The small nuclear effects in deuterium are neglected. 

In the EPS09 analysis, we have altogether 929 data points from 32 different data sets at our disposal. Three types 
of data are used: deep inelastic lepton-nucleus scattering (DIS), Drell-Yan (DY) dimuon production in p+A collisions 
from FNAL, and most lately - and included so far only in EPS08 (LO) E0) and in EPS09 (NLO, LO) Q - also 
neutral-pion production at mid-rapidity in d+Au collisions from BNL-RHIC. Note that unlike in EPS08, we do not 
utilize the BRAHMS forward-rapidity hadron data from d+Au in the EPS09 analysis. 

The goodness measure x 2 , which is minimized as described in Sec.[T| is here a generalized one 



X 2 ({a})^w NX %({a}), where x%({a}) ee ( 1 {^ a}) ) + £ 



N \ N ieN 



/ JV ({a})A-2H{o})-' 2 
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where we give an additional weight wn to those data sets N which provide important constraints but whose number 
of data points is so small that these constraints would otherwise escape unnoticed. This standard weighting procedure 
does not cause a bias as long as no significant tension between different data sets arises. The multiplier /jv({a}), 
which minimizes the x%({ a }) f° r a data set {D{\m at every round of iteration (needed only for the RHIC data now), 
is an output of the analysis, giving the best estimate of the overall normalization which is still consistent with the error 
p.norm pj-oyj^ed by the experiment. 

We minimize the global x 2 with respect to 15 fit parameters a. L , which quite obviously are correlated with each 
other. The minimum is found to be at Xq = 731.3 for 929 data points. Once Xq (the best fit) is determined, we perform 
an error analysis using the Hessian method. First, to find those linear combinations of ai which are uncorrected, we 
diagonalize the Hessian matrix H, whose elements are the second derivatives of x 2 at Xo- Then, denoting by Zi the 
suitably normalized parameter eigendirections, we can expand the x 2 as follows: 



X Xo ^ 2 d ai da, 

ij 



, = n 



ScuScLj + 0{{Sa t } 3 ) ^Xo + Y] Sa t H i:i 5a 3 = xl + V] z 2 , (4) 
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Figure 1: Elements in the EPS09 analysis. Left: A sketch of the EPS09 fit functions Rf(x) and the role of certain fit 
parameters. Right: \ 2 ~ Xo as a function of each linear combination zi, and its comparison with the ideal quadratic 



behaviour A\ 



(the shaded band). Bottom: The 90% confidence limits for each data set in the eigendirection 



z 10 . The filled boxes show the location of the \ 2 minimum for each data set, while the global minimum is at z = 0. 
All panels are from [ 1 1 . 



After this, the error propagation becomes more straightforward and we can obtain the error for a physical nPDF- 
dependent quantity X through the variations of the uncorrected parameters Zi as (for details, see £[]), 

where we have taken the same increase, Ax 2 = (Szi) 2 , in each eigendirection i of H. Figure [T| (right) shows that the 
quadratic approximation in Eq. M above is working well sufficiently near the \ 2 minimum. No unique procedure, 
however, exists for obtaining Ax^. We use the "90% confidence criterion" originally developed by CTEQ [22 J, and 
explained in detail in [1|. The idea, illustrated in Fig. [I] (bottom), is as follows: the \ 2 minimum for each data set of 
N data points is first determined separately for each eigendirection i. Then, we let the parameter Zi grow(decrease) 
from the extremum location up(down) to the point that the \ 2 of this data sets reaches the value below which the x 2 
should reside within a 90% probability, if an iV-dimensional Gaussian probability distribution for the \ 2 is assumed. 
Repeating this procedure for all the data sets, and in the end taking the stringest of all the obtained upper/lower limits, 
(Szf) 2 = Ax 2 {Szf) for each data set, and averaging over all eigendirections i, we arrive at 

, _ v A x 2 (Sz+) + A x 2 (5z7) (6z+) 2 + {Szjf 

= ^ 2N ^ 2N 5 °- (6) 

% i 

Whether A% 2 = 50 above reproduces the average size of the actual error bars in the data, needs of course be checked 
a posteriori. From Fig. [l] (right), however, we see that for this choice of A^ 2 , the quadratic approximation is still fine. 

Thus, the EPS09 nPDF package contains the central set (best fit) 5° and 15+15 error sets Sf, which are obtained 
by changing the fit parameters along the plus and minus directions of Zi so that the global x 2 grows by 50. Since the 
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error sets Sf 1 may shift a physical quantity X which the user wishes to study, into the same direction, we recommend 
the following prescription ll23ll for computing the upper and lower limits of X: 



(AX + ) 2 
(AX-) 2 



J2 [max {X(S+) - X(S ),X(Sr) - X(S ),0}] 2 , 

i 

£ [max {X(S°) X(Sf),X(S°) X(ST),0}] 2 ■ 



(7) 
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Figure 2: Left: The nuclear modifications Ry h , R^ h , i?g b and their uncertainties at Q 2 
From [ 1 1. Right: Q 2 evolution of these modifications at selected fixed values of x. 



1.69 GeV 2 and 100 GeV 2 . 



3 Results from EPS09 

Figure [2] (left) shows the average valence and sea quark and gluon modifications in a lead nucleus at the initial scale 
Qq and at a higher scale Q 2 = 10 GeV 2 according to the EPS09 NLO central set So (solid lines) and the error sets Sf 
(dotted lines). The shaded uncertainty band is computed using Eq. |7]i above. The largest uncertainties obviously reside 
at the smallest- x and largest- a; gluons, since these are the gluon regions worst constrained by the data. In the DGLAP 
evolution, Fig. |2](right), the small-a; gluon uncertainties, however, quickly shrink, while the large-x uncertainties not 
only remain but are also transferred into the large-x sea quarks. 

Figure [3] shows an example of the comparison of the EPS09 NLO results with the DIS data for the ratios 

* ^ ' ' 2 DIS/ ^ i—q,q,g 

(8) 

From this figure, and from a similar comparison of the DY cross section ratios (see Fig. 4 in Ref. [1 1), we confirm that 
with the choice A\ 2 = 50 the data uncertainties are transferred into the PDFs quite nicely: the error bands are indeed 
of the same size as the average error bars in the data. 

The effects of the DGLAP evolution in the nuclear modifications are perhaps best illustrated by Fig.|4]which shows 
the comparison of the EPS09 NLO results for F^/F^ (left) and for the DY cross section ratio (right) 



RU^M 2 )^ f^J ™ ~T > where <^ l r+X = E f!{M 2 )®ff(M 2 )®a^-+*{M 2 



(9) 
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Figure 3: Left: Rp 2 (x, Q 2 ) as computed in NLO using the EPS09 nPDFs (filled squares and error bands), and as 
measured by NMC 95 ll25l [261 (open symbols with error bars). Right: The same for R^ ls (x, Q 2 ) and SLAC data 
El. From H]. 



and x\ — (M 2 j \fs)e v , with M 2 being the invariant mass and y the rapidity of the lepton pair. The smallest-x panel 
on the left is currently the directmost constraint for nuclear gluons at x ~ 0(0.01). Smaller- a; DIS data at perturbative 
scales would be badly needed to further constrain the nuclear gluon distributions at smaller x. 

Comparison of the EPS09 NLO results for the inclusive mid-rapidity pion production ratio in minimum-bias d+Au 
and p+p collisions at RHIC, 
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i,j,k=q,q,g 

with the PHENIX data is shown in Fig. [5] (left). We use the KKP fragmentation functions [30] here, but we have 
checked that with the newer ones, AKK08 [31] or nDSS f32l[33l . the results would be the same. The NLO code by 
Aversa et al. [34| from the INCNLO ||351 compilation has been utilized after changing the order of integrations to speed 
up the computation (see [1]). Interestingly, and importantly, a good and essentially tensionless fit of these RHIC data 
is obtained simultaneously with the DIS and DY data. This data set, especially the downward trend at > 5 GeV, 
also offers a valuable constraint for the antishadowing and EMC effects in the gluon distributions. 

Figure|5]also compares the NLO results from the global DGLAP analyses HKN07 fl]D, nDS H3 and EPS 09. The 
main difference between these sets is the existence of the gluon antishadowing/EMC-effect structure present only in 
EPS09. This is mainly due to the fact that the PHENIX ir° data, which shows such behaviour, are included only in 
EPS09. The number of data points included in the EPS09/HKN07/nDS analyses is, correspondingly, 929/1241/420, 
while the obtained \ 2 per data point is, correspondingly, 0.79/1.2/0.72. 
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Figure 4: Left: The calculated NLO scale evolution (solid lines and error bands) of F 2 Sn /F 2 c and the NMC data [28 ] 
for fixed values of x. Right: R& Y (x u M 2 ) results from EPS09 NLO (filled squares and error bands) as a function of 
x\ compared with the E866 [27| data (open squares). From JT|. 
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Figure 5: Left: i?dAu from EPS09 NLO (solid line, blue error band) at y = for inclusive ir° production and the 
PHENIX 1 29 1 data (open squares). The error bars show the statistical and the yellow band the systematic errors. The 
prediction from HKN07 (NLO) is shown by the dashed line. The data have been multiplied by /jy- = 1.03, see Eq. <|3j, 
which is an output of the EPS09 analysis. The STAR data |36] (open circles), multiplied by — 0.90, are shown 
only for comparison. Right: The average valence and sea quark, and gluon modifications at Q 2 = 1.69 GeV 2 and 
Q 2 = 100 GeV 2 for Pb nucleus from the NLO HKN07 [18], nDS El and EPS09. From 0Q. 



4 EPS09 vs. BRAHMS forward-rapidity data 

As an example of using the EPS09 package, we consider the forward-rapidity hadron production, measurements for 
which have been done by BRAHMS ll37l . Figure [6] shows the invariant distributions EdN/d 3 p in p+p (left) and in 
d+Au collisions (right). The error bands shown in the p+p case have been obtained by using the central set together 
with the 40 error sets from CTEQ6.1M, and the ones for the d+Au case by using the central sets from CTEQ6.1M 
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Figure 6: Top left: EdN/d 3 p of negatively charged hadrons in p+p and d+Au collisions vs. px- The data (open 
squares) with error bars are from [ 37 1 with statistical and systematical errors added in quadrature. The error band is 
the free-proton PDF uncertainty given by CTEQ6.1M. Top right: The same for d+Au, except now the error bands 
are obtained using the error sets of both CTEQ6.1M and EPS09. Bottom: The corresponding ratio i?dAu- The error 
bands here correspond to the EPS09NLO error sets. From JT] 



and EPS09NLO together with their 40(CTEQ6.1M) + 30 (EPS09NLO) error sets as explained in Sec. 4 of Ref. (TJ. 
We notice that the data are systematically above the NLO pQCD calculation in the p+p case at 77 = 3.2, while in both 
nuclear cases considered, and also for 77 = 2.2 in p+p, there is an agreement. This discrepancy is why these BRAHMS 
data was excluded from the EPS09 analysis. After this, it is not a surprise that the data for the nuclear ratio i?dAu 
at T] = 3.2 and pt > 2 GeV, as shown by the bottom panel of Fig. [6] lie clearly below the EPS09NLO prediction. 
Thus, the NLO EPS09 global analysis suggests that the strong suppression seen in the BRAHMS i?dAu at 77 = 3.2 
and px > 2 GeV is due to an excess in the p+p data, not due to a suppression in the d+Au data. 



5 Summary and outlook 

The conclusion from the EPS09 NLO global DGLAP analysis is that there is an excellent agreement between NLO 
pQCD and the hard-process nuclear data for deeply inelastic lepton-nucleus scattering, DY-dilepton production in 
p+A, and central-rapidity tt° production in d+Au collisions in the nPDFs' kinematical range 0.005 < x < 1, 1.69 < 
Q 2 ^ 150 GeV 2 . The goodness of the global fit with 929 data points is x 2 /N = 0.79, and remarkably, no significant 
tension between the data sets from different types of processes is found. Collinear factorization seems to work well in 
describing high-energy inclusive nuclear hard processes. 

The EPS09 nPDF package contains the central set (best fit) + 30 error sets, both in NLO and LO, and they are 
downloadable at https://www.jyu.fi/fysiikka/en/research/highenergy/urhic/nPDFs. With these error sets, it is now, 
finally, possible for anybody to study how the nPDF-uncertainties propagate into nuclear hard cross sections. This 
should be particularly useful for the precision analysis of hard probes as QCD matter signatures, and also, e.g., for 
planning detector upgrades in the near future. 

Further tests of factorization and nPDFs will be provided e.g. by RHIC data for direct photons in d+Au, and 
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both RHIC and LHC data in Au+Au collisions. The corresponding pQCD cross sections have been discussed e.g. in 
Ref. £38), and also EPS09NLO-related direct photon studies both for RHIC and LHC have been launched (391. To 
constrain the nuclear quark distributions better, also (anti)neutrino-nucleus DIS data should be added into the global 
analysis of the nPDFs. Recent studies [40 1, suggest that the agreement with the CHORUS, CDSHSW and also most 
of the NuTeV data (see [40] for references) is in fact quite good. Before including the neutrino-DIS data into the 
EPS09-type global analysis consistently, however, one needs to extend the analysis to a general-mass framework. 

Further tests of factorization will also be given by heavy-quark production and forward-rapidity pion data in d+Au 
collisions at RHIC. To properly test the applicability of factorization at the available highest energies as well - and to 
get the most reliable comparison baselines for the nuclear hard processes - it would quite obviously be very important 
to have the p+Pb runs scheduled at the LHC, too. For recent studies of constraining the nPDFs at the LHC, see e.g. 
Refs. I4T1I421 . The most direct constraints for the still unconstrained lowest- a; nuclear gluons one would, however, get 
from the planned future colliders eRHIC and in particular LHeC. 

Finally, regarding a longer-term outlook for the global analyses, in order to obtain a minimal number of uncorre- 
lated error sets simultaneously for the free-proton PDFs and for the nPDFs, one should eventually combine the global 
analyses for the free and bound proton PDFs into one single master analysis. 

Acknowledgements: KJE thanks the Vilho, Yrjo and Kalle Vaisala fund of the Finnish Academy of Science and 
Letters, and the Academy of Finland (Project nr. 133005) for financial support. 
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